#ifndef __ALLCODE_PHYSIKS_H_
#define __ALLCODE_PHYSIKS_H_


#include <vector>
#include <iostream>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_const_mksa.h>
#include <GL/glut.h>


using namespace std;


typedef unsigned int    uint;
typedef double  PHYSIKSCALAR;



namespace PhysiksBase {

/** \class Vector
 * A simple wrapper class around gsl/gsl_vector. Without this, one would have an unpleasant
 * blend of C-style code alongside C++ code. Note that it is work in progress, with new
 * functions added when they are needed. */
class Vector {
    public:
        /** \brief  Constructor 1.
         * \param size The dimension of the vector
         * \param init the initialisation value of all vector elements */
        Vector(unsigned int size, double init = 0.0);

        /** \brief  Constructor 2.
         * \param init the initialisation gsl_vector */
        Vector(const gsl_vector *init);

        /** \brief  Copy Constructor.
         * \param init the reference Vector */
        Vector(const Vector &init);

        virtual ~Vector();


        /** \brief  In place scalar multiplication (i.e. with side effect).
         * \param fac the scaling factor */
        void scale(double fac);

        /** \brief  Set the vector to the zero vector. */
        void set_zero();

        /** \brief  Set the value of a coordinate. No bounds checking.
         * \param index the index of the coordinate (0-based)
         * \param value the value to change the coordinate to */
        void set(unsigned int index, double value);

        /** \brief Sets the coordinates of this vector to those of the parameter.
         * \param copy the template vector, of which coordinates are to be copied over  */
        void set(const Vector &copy);
        
        /** \brief  The Euclidean norm (magnitude) of the Vector.
         * \return the euclidean norm. */
        double norm() const;

        /** \brief  A unit vector in the same direction as the calling vector.
         * \return a unit vector in the same direction as the calling vector */
        Vector normalize() const;
        
        /** \brief Effectively "clones" the current vector.
         * \return  an identical copy of this vector    */
        Vector duplicate() const;

        /** \brief  Returns the dimension of the vector.
         * \return calling vector's dimension */
        unsigned int dim() const;

        /** \brief  Get accessor for an arbitrary coordinate. Bounds checking offered.
         * \param index The coordinate's index (0-based).
         * \return the value of the respective coordinate. 0 if index is out of bounds */
        double get(unsigned int index) const;
        
        /** \brief  Vector cross product. Does not check dimensions of cross-product operands, so
         *          vectors of dimension other than 3 might yield non-sensical results.
         * \param   rhs the 2nd operand of the cross product.
         * \return  the resultant vector, orthogonal to both the caller and rhs. */
        Vector cross(const Vector &rhs) const;


        /** Vector addition overload. */
        friend Vector operator+(const Vector &lhs, const Vector &rhs);
        /** Vector addition & assignment overload. */
        friend void operator+=(Vector &lhs, const Vector &rhs);
        /** Vector subtraction & assignment overload. */
        friend void operator-=(Vector &lhs, const Vector &rhs);
        /** Vector subtraction overload. */
        friend Vector operator-(const Vector &lhs, const Vector &rhs);
        /** Vector negation overload. */
        friend Vector operator-(const Vector &rhs);
        /** Vector standard inner product (dot product) overload. */
        friend double operator*(const Vector &lhs, const Vector &rhs);
        /** Vector scalar multiplication overload. */
        friend Vector operator*(const Vector &lhs, double fac);
        /** Vector scalar multiplication overload. */
        friend Vector operator*(double fac, const Vector &rhs);
        /** Vector indexing (random coordinate access) function. No bounds checking. */
        double operator[](unsigned int k) const;

        /** \brief  Returns a pointer to the gsl_vector core of the calling Vector.
         * PLEASE DO NOT USE THIS. It exists solely for (rare) internal usage. */
        gsl_vector *core() const;

    private:
        /** The dimension of the Vector. This should be immutable, no? */
        unsigned int dimension;
        
        /** The gsl_vector core. */
        gsl_vector *cont;
};

}   // End namespace PhysiksBase.




namespace PhysiksBase {

class NewtonianGravitationalField;
class CoulombAccelerationField;

/** \class PhysicalObject 
 * \brief The fundamental dynamic element of a simulation. */
class PhysicalObject {
    public:
        double mass;
        double charge;
        Vector position;
        Vector velocity;
        Vector net_force;
        
        /** Flag specifying whether or not this object should be fixed at a point in Rn */
        bool immovable;
        
        /** Gravitational field created by this object  */
        NewtonianGravitationalField gravitational_field;
        
        /** Coulomb (electric) acceleration field created by this object   */
        CoulombAccelerationField coulomb_acceleration_field;

        /** \brief  Default parameterless constructor. Mass and charge are initialised to 0, and
         * all vectors are initialised to the zero vector. */
        PhysicalObject();
        
        /** \brief  Explicit constructor.
         * \param   mass    the object's mass.
         * \param   charge  the object's charge.
         * \param   position    the object's initial position in n-space.
         * \param   velocity    the object's initial velocity in n-space. */
        PhysicalObject(double mass, double charge, const Vector &position,
                       const Vector &velocity);
        
        /** \brief Draws the OpenGL representation of this type of PhysicalObject. */
        virtual void draw() const;
        
        virtual ~PhysicalObject();
};

//=============================================================================

class ExtendedObject : public PhysicalObject {
    public:
        Vector angular_velocity;
        
        /** \brief  Default parameterless constructor. In addition to the initialisations of 
         * PhysicalObject's parameterless constructor, angular_velocity is initialised to the zero.
         * vector. */
        ExtendedObject();
        
        /** \brief  Explicit constructor.
         * \param   angular_velocity    the extended object's angular velocity. */
        ExtendedObject(double mass, double charge, const Vector &position, const Vector &velocity);
        
        virtual ~ExtendedObject();
};

//=============================================================================

class Sphere : public ExtendedObject {
    protected:
        double radius;
        
        /** This sphere's constant surface elasticity. Used primarily for collision reaction */
        double elasticity;
        
    public:
        /** \brief  Default parameterless constructor. In addition to the initialisations of
         * ExtendedObject's parameterless constructor, radius is initialised to 1. */
        Sphere();
        
        /** \brief  Explicit constructor.
         * \param   radius  the spherical object's radius. */
        Sphere(double mass, double charge, const Vector &position, const Vector &velocity, 
               double radius);
        
        virtual ~Sphere();
        
        /** \brief Draws a sphere at the appropriate location in GL-space, representing this
         *         instance. */
        void draw() const;
        
        double get_radius() const;
        void set_radius(double radius);
        double get_elasticity() const;
        void set_elasticity(double val);
};

//=============================================================================

class Plane : public ExtendedObject {
    protected:
        /** \brief  Flag controlling the nature of the plane's interactions with other objects.
         * 
         * If true, the plane is taken to be infinite (i.e. the set of ALL linear combinations of
         * the two basis vectors. Otherwise, it is restricted to the set of all linear combinations
         * of the two basis vectors with coefficients <= 1. In this way, the length and breadth of
         * the bounded plane are specified by the magnitudes of basis_vec1 and basis_vec2. */
        bool unbounded;

    public:
        /** The plane's unit normal vector. */
        Vector *normal_vec;

        /** One of the plane's two basis vectors. */
        Vector *basis_vec1;
        
        /** The second of the plane's two basis vectors. */
        Vector *basis_vec2;
        

        Plane(const Vector &normal, const Vector &origin);

        /** \brief  Explicit constructor.
         * \param   vec1    Basis vector 1 for the plane.
         * \param   vec2    Basis vector 2 for the plane.
         * \param   origin  the common base of vec1 and vec2.   */
        Plane(const Vector &vec1, const Vector &vec2, const Vector &origin);
        
        virtual ~Plane();
        
        void draw() const;
        
        /** \brief  Set the plane's unbounded flag.
         * \param   flag    true if the plane must be bounded, false otherwise. */
        void bounded(bool flag);
        
        /** \brief  Returns the value of the plane's ``bounded'' flag
         * \return  True if the plane is bounded, false otherwise. */
        bool bounded() const;
        
        bool contains(const Vector &vec) const;
        
        
        /** Returns a unit vector normal to this plane  */
        Vector normal() const;
        
        /** Sets this plane's normal vector */
        void set_normal(const Vector &norm);
};



}   // End of namespace PhysiksBase



namespace PhysiksBase {


/** A base class for all fields in Physiks. */
class Field {};

/** A common interface for all vector fields. */
class VectorField : public Field {
    public:
        /** \brief Returns the corresponding output vector, given a position in Rn */
        virtual Vector operator()(const Vector &) const = 0;
};

/** A common interface for all scalar fields. */
class ScalarField : public Field {
    public:
        /** \brief Returns the corresponding output scalar, given a position in Rn */
        virtual double operator()(const Vector &) const = 0;
};


/** Built-in magnetic (B) field.    */
class MagneticField : public VectorField {};

/** Built-in electric (E) field.    */
class ElectricField : public VectorField {};

/** Built-in force (F) field.   */
class ForceField : public VectorField {
    public:
        virtual Vector operator()(const Vector &) const = 0;
};

/** Built-in acceleration field.    */
class AccelerationField : public VectorField {};

/** Built-in velocity field.    */
class VelocityField : public VectorField {};


/** Built-in representation of a Newtonian gravitational acceleration field.   */
class NewtonianGravitationalField : public AccelerationField {
    private:
        
        /** \brief The center of the field, from/to which force vectors emanate.   */
        PhysicalObject *center;
        
    public:
        
        /** \brief Constructor accepting the object to be bound to as the sole parameter.
         * \param   obj the object at the center of the field, responsible for the attractive
         *              force.  */
        NewtonianGravitationalField(PhysicalObject &obj);
        
        /** \brief  Implementation of the acceleration field function.
         * \param   pos a position in Rn
         * \return  the <b>acceleration</b> vector corresponding to the given position  */
        Vector operator()(const Vector &pos) const;
        
        virtual ~NewtonianGravitationalField();
};

/** Built-in representation of a Coulomb electric acceleration field.   */
class CoulombAccelerationField : public AccelerationField {
    private:
        
        /** \brief The center of the field, from/to which force vectors emanate.   */
        PhysicalObject *center;
        
    public:
        
        
        /** \brief Constructor accepting the object to be bound to as the sole parameter.
         * \param   obj the object at the center of the field, responsible for the attractive
         *              force.  */
        CoulombAccelerationField(PhysicalObject &obj);
        
        /** \brief  Implementation of the acceleration field function.
         * \param   pos a position in Rn
         * \return  the <b>acceleration</b> vector corresponding to the given position  */
        Vector operator()(const Vector &pos) const;
        
        virtual ~CoulombAccelerationField();
};

}   // End of namespace PhysiksBase




namespace PhysiksBase {


/** A World encapsulates all information about fields, objects and parameters of the physical
 * environment. Force fields and objects (among other things) are registered in a specific
 * world instance, and integration is initiated via the instance itself. */
class World {
    public:
        /** the dimension of all tuples used for this world. */
        const int dimension;
        
        /** time step used for integration. Not constant to enable runtime modification. */
        double timestep;
        
        /** flag controlling whether or not data should be output (currently only to stdout) */
        bool output;

        /** \brief  Constructor 1.
         * \param dimension the dimension of all tuples used for this world
         * \param timestep the initial timestep for integration */
        World(int dimension, double timestep);

        /** \brief  Constructor 2.
         * \param init_objects An initialisation reference for this world's object table. */
        World(int dimension, double timestep, vector<PhysicalObject *> &init_objects);

        /** \brief  Constructor 3.
         * \param   init_force_fields   An initialisation reference for this world's force field
         *                              table. */
        World(int dimension, double timestep, vector<PhysicalObject *> &init_objects,
              vector<ForceField *> &init_force_fields);
        
        virtual ~World();

        /** \brief  Adds an object to the calling world's list of active objects.
         * \param obj the object to be added */
        void add_object(PhysiksBase::PhysicalObject &obj);
        
        void add_single_object_proc(void (*)(PhysicalObject *));
        void add_double_object_proc(void (*)(PhysicalObject *, PhysicalObject *));


        /** \brief  Commence integration for the calling world.
         * \param start_time the effective absolute time in the world when integration begins
         * \param end_time the effective absolute time in the world when integration stops */
        void integrate(double start_time, double end_time);
        
        /** \brief Performs a single iteration of the timestep and updates all world parameters. */
        void step();

        /** \brief Returns the current absolute time in the world.
         * \return the current absolute time */
        double current_time() const;
        
        /** \brief References the object at the specific index in the object table.
         * \param   index   the index of the queried object in this world's object table
         * \return  a pointer to the queried object. <b>The object is modifiable through the
         *          returned pointer.</b>    */
        PhysicalObject *object(unsigned int index) const;
        
        /** \brief  Returns the number of objects present in this world.
         * \return  the number of active objects in the world.  */
        unsigned int num_objects() const;
        
        // This is one way of implementing custom force fields, but it forces the user to muck
        // around with the internals of the argument. Not good.

        // User-defined unary procedures.
        vector<void (*)(PhysicalObject *)> single_object_procs;
        // User-defined binary procedures.
        vector<void (*)(PhysicalObject *, PhysicalObject *)> double_object_procs;

    private:
        /** A table of all objects in the world. */
        vector<PhysicalObject *> objects;

        /** A table of all defined force fields in the world. */
        vector<ForceField *> force_fields;
        
        /** A table of all defined accleration fields in the world. */
        vector<AccelerationField *> acceleration_fields;
        
        /** A table of all defined velocity fields in the world. */
        vector<VelocityField *> velocity_fields;
        
        /** \brief The current absolute time in the world. Will need to be discarded if someone
         *         decides to add relativistic computations to PhysiKs' capabilities. */
        double time;
        
        /** \brief Applies the user defined custom procedures for this world. */
        void apply_user_defined_procs();

        /** \brief Applies the user defined fields to all objects in the world's object table. */
        void apply_user_defined_fields();
        
        /** Please ignore. Needs to be deleted. */
        void print_data();
};


}   // End of namespace PhysiksBase




namespace PhysiksBase {
    namespace Computation {
        /** \brief Applies each force field specified in the 2nd parameter to all objects in the
         *         given object table.
         * \param objTable   an object table containing pointers to the objects to be updated
         * \param force_fields   all force fields to be applied to each object in objTable  */
        void apply_force_fields(vector<PhysicalObject *> &objTable,
                                const vector<ForceField *> &force_fields);
        
        /** \brief Applies each acceleration field specified in the 2nd parameter to all objects in
         *         the given object table.
         * \param objTable   an object table containing pointers to the objects to be updated
         * \param force_fields   all acceleration fields to be applied to each object in objTable
         */
        void apply_acceleration_fields(vector<PhysicalObject *> &objTable, 
                                       const vector<AccelerationField *> &accel_fields);

        /** \brief Applies each velocity field specified in the 2nd parameter to all objects in the
         *         given object table. */
        void apply_velocity_fields(vector<PhysicalObject *> &objTable, 
                                   const vector<VelocityField *> &vel_fields);
        
        /** \brief Computes the distance between two physical objects in Cartesian-Euclidean space.
         * @param a object a
         * @param b object b
         * @return the scalar distance between a and b */
        double compute_distance(const PhysicalObject &a, const PhysicalObject &b);
        
        Vector compute_connecting_vector(const PhysicalObject &a, const PhysicalObject &b);

        /** \brief Sets the force vectors of all objects in objTable to the zero vector.
         * @param objTable the object table containing pointers to the relevant objects */
        void set_force_vectors_zero(vector<PhysicalObject *> &objTable);

        /** \brief Updates the force vectors of all objects in objTable based on built-in force
         *         fields, and user-defined fields. */
        void update_forces(vector<PhysicalObject *> &objTable);

        /** \brief Updates the velocity vectors of all objects in objTable based on their force
         *         vectors.
         * @param time_step the time interval in the equation v = v_0 + a*t */
        void update_velocities(vector<PhysicalObject *> &objTable, double time_step);

        /** \brief Updates the position vectors of all objects in objTable based on their net force
         *         and velocity vectors.
         * @param time_step the time interval in the equation s = s_0 + v_0*t + 0.5*a*t^2 */
        void update_positions(vector<PhysicalObject *> &objTable, double time_step);

        /** \brief Computes the Newtonian gravitational force between two objects.
         * @return the magnitude of the gravitational force between obj1 and obj2 */
        double compute_gravitational_force(const PhysicalObject &obj1, const PhysicalObject &obj2);

        /** \brief Calculates the signed magnitude of the coulomb force between two charges.
         * @param a Object a
         * @param b Object b
         * @return the signed magnitude of the Coulomb force between charges a & b */
        double compute_coulomb_force(const PhysicalObject &obj1, const PhysicalObject &obj2);
    }
}



namespace PhysiksBase {

namespace Collisions {
    
    /** Detect collisions among various PhysicalObjects within a world. */
    bool detect_collisions(const World &world);
    
    /** Check for a collision between the given point and sphere. */
    bool detect_point_sphere(const PhysicalObject &obj, const Sphere &sph);
    
    /** Check for a collision between the given point and plane. */
    bool detect_point_plane(const PhysicalObject &obj, const Plane &plane);
    
    /** Check for a collision between two spheres. */
    bool detect_sphere_sphere(const Sphere &sph1, const Sphere &sph2);
    
    /** Check for a collision between the given sphere and plane. */
    bool detect_sphere_plane(const Sphere &sphere, const Plane &plane);
    
    
    /** Each of the following functions adjusts the force vectors of its
     * parameters based on the relevant collision information (e.x. positions
     * of objects relative to each other, elasticity of spheres, etc.)  */
    void react_point_sphere(PhysicalObject &obj, Sphere &sph);
    void react_point_plane(PhysicalObject &obj, Plane &plane);
    void react_sphere_sphere(Sphere &sph1, Sphere &sph2);
    void react_sphere_plane(Sphere &sphere, Plane &plane);
    
}   // End of namespace Collisions
}   // End of namespace PhysiksBase



#endif
